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ABSTRACT 

Satellite  ocean  color  remote  sensing  is  plagued  by  loss  of  coverage  due  to  cloud  obscuring,  glint  contamination,  at¬ 
mospheric  correction  failures,  and  other  issues.  We  have  developed  a  simple  and  efficient  technique  for  estimating 
missing  remote  sensing  data  by  taking  advantage  of  the  inter-pixel  spatial  and  temporal  coherency  of  individual 
ocean  color  products.  The  technique  first  employs  a  limited  iterative  triangular  interpolation  procedure.  This 
procedure  attempts  to  select  three  neighboring  pixels  forming  the  tightest  triangle  enclosing  the  data  point  we 
are  attempting  to  recover;  and  then  interpolating.  On  failure  to  find  three  suitable  neighbors,  a  second  procedure 
is  employed  which  attempts  to  recover  missing  data  points  by  using  a  time  dependent  “latest  pixel”  replacement. 
This  procedure  replaces  the  missing  data  point  with  the  most  recent  data  point  collected  at  that  grid  point  within 
the  last  seven  days.  This  technique  has  been  applied  to  MODIS  (MODerate  resolution  Imaging  Spectrometer) 
ocean  color  products  of  phytoplankton  absorption,  back-scattering  coefficient,  and  chlorophyll  concentration  to 
produce  cloud  free  bio-optical  products  on  a  daily  basis  and  provide  a  new  capability  for  monitoring  coastal 
processes.  We  demonstrate  a  new  method  on  MODIS  products  and  show  how  bio-optical  properties  change  over 
a  daily  and  monthly  time  scale. 

Keywords:  Coastal,  Ocean  Color,  Remote  Sensing,  Composite,  Cloud  Removal 

1.  INTRODUCTION 

Satellite  remote  sensing  provides  the  benefits  of  wide  spatial  coverage  and  high  temporal  revisitation.  This 
provides  an  advantage  over  other  remote  sensing  platforms  such  as  aircraft.  The  nature  (it’s  in  space)  and 
in-flexibility  (fixed  look  angle  and  time  of  over  flight)  of  satellite  remote  sensing  platforms  cause  them  to  be 
vulnerable  to  a  few  issues  that  other  platforms  can  more  easily  avoid;  namely  cloud  obscuring  and  sun  glint. 
These  two  issues  along  with  other  common  causes  of  data  coverage  loss  such  as  atmospheric  correction  failure  or 
product  algorithm  failure  are  a  constant  plague  to  the  field  of  remote  sensing. 

The  open  and  coastal  ocean  processes  occur  at  different  time  and  space  scales.  Satellite  ocean  color  imagery 
captures  the  near  surface  properties  of  these  processes  at  specific  “windows”  of  opportunities  based  on  the  revisit 
time  of  the  satellite  and  ability  to  see  the  ocean  surface  (cloud  free,  glint  free,  etc.).  The  revisit  time  for  MODIS 
is  fixed,  however  attempts  to  characterize  the  ocean  processes  are  often  limited  as  a  result  of  clouds,  producing 
data  voids.  These  data  voids  hinder  our  ability  to  interpret  the  ocean  processes.  Specifically,  the  dynamic 
movements  of  ocean  features  are  disrupted  by  these  data  voids  and  it  is  difficult  to  link  the  bio-optical  properties 
as  a  response  to  circulation.  Typically,  to  understand  the  changing  structure,  a  7-day  latest  pixel  composite,1 
is  used  to  fill  in  the  data  voids.  Here  the  most  recent  valid  pixel  is  used  to  represent  the  ocean  surface.  If  the 
pixel  is  not  valid,  the  previous  day’s  value  is  used.  If  that  day  is  not  valid  then  the  previous  day’s  pixel  is  used. 
This  procedure  goes  back  for  one  week  trying  to  find  a  valid  pixel.  If  no  valid  pixel  is  found  for  a  week,  the 
value  is  not  determined  and  the  area  remains  black.  Figure  1  shows  an  example  of  the  Northern  Gulf  of  Mexico 
indicating  no  retrieved  data  as  black  areas.  Using  this  procedure,  a  daily  image  can  be  produced  to  represent 
the  most  recent  satellite  observations  of  the  bio-optical  properties.  This  procedure  works  well  in  open  mesoscale 
processes  which  do  not  change  too  rapidly  and  where  ocean  processes  have  a  residence  time  of  ~  3  —  4  days  for 
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Figure  1.  The  latest  pixel  composite  (LPC)  for  the  Northern  Gulf  of  Mexico  for  MODIS  Aqua  surface  chlorophyll  for  July 
18  with  previous  compositing  beginning  July  12.  The  black  areas  indicate  that  for  the  last  seven  days  no  valid  pixels  were 
available. 


observable  change.  The  dynamics  of  ocean  processes  become  less  stable  in  coastal  areas  as  the  processes  become 
more  dynamic  with  resident  times  approaching  hours  (tides,  winds,  etc.).  For  these  coastal  areas  a  more  robust 
method  is  required  to  better  represent  the  optimum  replacement  of  “data  voids” . 

To  improve  the  satellite  chlorophyll  product,  we  have  developed  a  (fairly)  simple  and  (somewhat  computation¬ 
ally)  efficient  technique  for  estimating  missing  remote  sensing  data  points  by  taking  advantage  of  the  inter-pixel 
spatial  and  temporal  coherency  of  individual  ocean  color  products.  We  make  the  assumption  that  ocean  color 
products,  on  any  grid  scale,  possess  a  value  similar  in  magnitude  to  their  neighbors  in  nearby  grid  locations.  We 
make  a  similar  assumption  in  temporal  space.  These  assumptions  allow  us  to  use  these  nearby  data  points  to 
regenerate  our  missing  data  points. 


2.  MODIS  PROCESSING 

MODIS  is  a  multi-spectral  sensor  with  36  discrete  spectral  bands.  MODIS  has  nine  bands  at  1km  resolution  cov¬ 
ering  the  visible  spectrum  between  412nm  and  869nm,  three  bands  at  500m  and  two  bands  at  250m.  The  MODIS 
data  was  processed  using  the  Automated  Processing  System2  Version  3.6  developed  by  the  Naval  Research  Labo¬ 
ratory  (NRL).  The  Automated  Processing  System  (APS)  is  a  collection  of  software  programs  assembled  by  NRL 
which  handles  the  complete  processing  of  satellite  imagery  after  reception  including  atmospheric  correction,  the 
application  of  bio-optical  algorithms  to  derive  optical  properties,  and  warping  of  the  imagery  to  a  standard  map. 
APS  implements  a  modification  of  the  NASA  SeaDAS  (v5.1)  atmospheric  correction  procedure.3’4 

APS  software  identifies  non- valid  ocean  pixels  by  screening  for  1)  Clouds  2)  Glint  contamination  3)  negative 
radiances  from  errors  in  atmospheric  correction  i.e.  epsilon  failures  4)  bio-optical  algorithm  failures.  Data  voids 
are  numerous  and  unavoidable  and  extremely  problematic  in  certain  ocean  regions.  These  voids  can  limit  the 
utility  of  the  satellite  data  especially  for  coastal  applications. 

APS  implements  a  number  of  traditional  compositing  techniques  such  as  weekly,  monthly,  and  yearly  averages. 
Pixels  which  have  been  identified  as  voids  are  excluded  from  the  compositing  methods.  Additionally,  APS 
implements  a  “latest  pixel”  composite  (LPC).  This  composite  method  attempts  to  replace  missing  data  points 
with  data  points  at  the  same  grid  location  collected  in  the  recent  past.  This  type  of  composite  preserves  good 
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Figure  2.  Example  of  how  the  triangular  interpolation  procedure  works.  The  data  voids  (Black)  sire  filled  in  using  radius 
distances  (referred  to  as  layers  a,  b,  c,  ...). 

data  points,  since  a  good  data  point  is  not  mixed  (averaged)  with  any  other  data  point.  The  “7-day  latest  pixel” 
composite  will  search  up  to  seven  days  into  the  past  to  find  a  good  data  point. 

3.  METHODS  AND  RESULTS 

The  triangular  compositing  technique  initially  employs  a  limited  iterative  triangular  interpolation  procedure. 
This  procedure  attempts  to  select  three  neighboring  pixels  forming  the  tightest  triangle  enclosing  the  data  point 
we  are  attempting  to  regenerate.  The  radius  of  the  search,  out  from  the  center  of  the  missing  data  point,  is 
restricted  to  ten  grid  points  (or  ten  layers,  see  figure  2).  Once  the  three  data  points  are  found,  a  triangular 
interpolation  is  performed,  and  the  result  becomes  the  new  value  for  the  missing  data  point.  If  the  three  data 
points  can  not  be  found,  then  the  value  is  left  unchanged.  Care  is  taken  to  avoid  interpolating  across  obstructions 
such  as  land. 

Figure  2  represents  an  example  of  the  search  procedure.  Dark  colored  circles  represent  data  voids,  and  light 
colored  circles  represent  data  points  with  known  data  values.  The  search  is  performed  using  a  layered  approach. 
For  the  case  of  the  dark  circle  marked  “1”,  the  first  layer  marked  “a”,  is  made  up  of  the  eight  circles  immediately 
surrounding  circle  “1”.  These  eight  data  points  are  added  to  a  list  of  data  points  and  for  every  combination  of 
three  data  points,  a  triangle  is  created  and  a  test  is  conducted  to  determine  whether  data  point  “1”  falls  inside 
the  triangle.  If  the  search  fails  at  layer  “a”,  then  the  search  is  expanded  to  include  the  data  points  from  the  next 
layer  out:  layer  “b”.  The  data  points  from  layer  “b”  are  added  to  the  list  of  data  points  which  already  include 
those  from  layer  “a”,  and  again  a  triangle  enclosing  our  target  data  point  is  searched  for.  This  procedure  is 
continued  until  three  valid  points  forming  a  triangle  around  our  point  are  found  or  the  maximum  search  radius 
is  reached.  The  red  circles  labeled  “1”  represent  the  three  points  selected  for  interpolation  forming  a  triangle 
around  point  “1”.  The  blue  circles  labeled  “2”  represent  the  three  points  selected  for  interpolation  forming  a 
triangle  around  point  “2” . 

We  have  the  ability  to  iterate  so  that  the  triangular  interpolation  procedure  is  performed  repeatedly.  Each 
successive  iteration  is  able  to  make  use  of  previously  interpolated  points  when  searching  for  the  three  neighboring 
pixels  for  the  current  triangular  interpolation.  The  iteration  can  be  performed  a  fixed  number  of  times  or 
repeatedly  until  no  additional  pixels  can  be  interpolated.  The  iteration  method  can  be  used  to  gradually  rebuild 
a  section  of  an  image  and  assists  in  reducing  the  neighbor  search  time  by  producing  similar  results  with  a  smaller 
search  radius. 
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Figure  3.  The  original  MODIS  Aqua  -  surface  Chlorohpyll  image  illustrating  the  data  voids  in  black. 

A  nearest  neighbor  search  is  also  implemented  els  a  fall  back  method.  As  its  name  suggests  this  is  simply  a 
search,  limited  by  the  search  radius,  for  the  neighboring  data  point  nearest  to  the  data  point  we  are  attempting 
to  regenerate. 

Additionally,  we  make  use  of  the  image  compositing  capabilities  of  APS.  We  use  APS’s  software  to  create 
the  “7-day  latest  pixel”  composite  described  above.  We  also  make  use  of  a  feature  which  allows  specifying  a 
limit  on  the  number  of  pixels  that  are  averaged  together  to  produce  a  data  point.  For  example,  if  five  scenes 
were  supplied,  and  a  valid  non-obscured  data  point  existed  in  each  scene,  then  five  pixels  would  be  averaged 
together  to  produce  the  composited  data  point.  If  we  limit  the  composite  to  two  pixels,  then  only  the  first  two 
encountered  data  points  will  be  averaged,  the  others  would  be  ignored. 

The  original  MODIS  Aqua  image  for  the  Northern  Gulf  of  Mexico  is  shown  in  figure  3.  The  black  areas 
offshore  are  clearly  flagged  as  clouds  and  the  coastal  areas  have  failed  due  to  atmospheric  correction  errors  which 
produced  negative  radiances.  This  image  is  contrasted  to  the  LPC  image  (figure  1)  where  a  small  amount  of  the 
black  areas  have  been  filled  in.  In  this  example  the  LPC  is  not  the  optimum  since  the  previous  seven  days  have 
collected  few  valid  pixels  in  the  region. 

The  new  methods  for  filling  in  the  void  pixels  on  this  image  are  performed  in  five  steps. 

1.  The  initial  stage  is  to  apply  a  triangular  interpolation  method  to  each  data  void,  searching  the  necessary 
layers  and  iterating  ten  times.  This  is  perhaps  the  most  important  step  where  the  spatial  interpolation 
fills  in  the  majority  of  the  data  drop  outs.  The  procedure  propagates  through  the  sequential  data  points 
filling  in  the  data  voids  using  a  search  radius  of  ten  pixels  or  layers.  Some  of  the  data  voids  will  remain 
after  this  procedure  is  performed  on  the  entire  image.  The  procedure  is  repeated  ten  times  to  continue  to 
replace  the  data  voids.  The  result  shown  in  figure  4  clearly  removes  the  cloud  noise  but  still  has  large  gaps 
in  coastal  areas. 

2.  The  next  step  is  to  use  a  nearest  neighbor  search  with  search  radius  one  to  create  a  border  around  known 
pixels  (figure  5).  This  procedure  is  similar  to  a  dilation  of  known  data  values,  and  enlarges  and  extends 
known  values  into  data  void  regions  for  bordering  pixels.  This  is  especially  useful  in  near  shore  coastal 
regions  and  forms  a  transition  of  the  bordering  data  with  the  next  step. 
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Figure  4.  Step  1  of  triangular  interpolation  of  the  data  voids  shows  replacement  of  the  majority  of  the  cloud  noise  and 
data  drop  outs. 
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Figure  5.  Step  2  creates  a  border  around  known  data  values  to  replace  data  voids;  especially  useful  as  a  transition  value 
when  combined  with  the  LPC  in  the  next  step. 
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Figure  6.  Composite  of  the  7-day  LPC  with  the  triangular  interpolation  methods  (steps  1  and  2). 


3.  The  next  procedure  is  to  create  an  “average”  composite  with  a  two  pixel  limit,  of  the  results  of  steps  1  and 
2  and  the  7-day  latest  pixel  composite.  This  will  produce  an  image  (figure  6)  which  will  retain  the  results 
produced  by  step  1,  and  which  will  average  the  results  of  step  2  and  the  7-day  latest-pixel  composite  where 
they  overlap.  The  7-day  latest-pixel  composite  will  fill  in  areas  for  which  there  is  no  overlap. 

4.  The  next  step  again  applies  the  triangular  interpolation  method  discussed  in  step  1  on  the  previous  step  3. 
This  image  has  now  included  additional  temporal  data  voids  from  the  previous  seven  days.  The  procedure 
again  iterates  ten  times,  similarly  to  Step  1.  The  result  (figure  7)  shows  increased  data  inclusion  in  the 
West  and  in  the  Lake  Pontchartrain.  Note  that  the  7-day  LPC  did  not  show  substantial  increased  data  as 
is  shown  in  figure  1. 

5.  The  final  step  is  used  to  remove  all  data  voids  and  completely  fill  the  data  with  the  best  estimate  of  the 
chlorophyll.  (In  certain  instances  a  data  value  is  required,  for  example  in  modeling  the  advection  field 
or  coastal  monitoring  where  our  best  guess  of  the  value  is  necessary).  An  estimate  of  the  data  points  at 
every  remaining  grid  point  is  produced  by  using  a  nearest  neighbor  search  with  search  radius  one  pixel, 
and  iterating  until  the  image  is  completely  filled  in  (figure  8).  These  continuously  replicated  values  can 
produce  significant  errors.  Note  that  the  data  replications  in  the  western  portion  of  the  image  have  high 
uncertainty.  Similarly  the  western  side  of  Lake  Pontchartrain  has  pixel  replication. 

4.  CONCLUSION 

A  MODIS  scene  containing  cloud  obscured  pixels  for  the  Mississippi  Bight  region  from  MODIS  Aqua  on  July 
18,  2007,  was  used  to  illustrate  the  procedure.  Here  we  followed  the  above  enumerated  procedure  to  replace  the 
data  voids  for  the  Chlorophyll  Concentration. 

We  have  applied  our  compositing  technique  to  additional  MODIS  ocean  color  products  of  phytoplankton 
absorption,  back-scattering  coefficient,  and  chlorophyll  concentration  to  produce  cloud  free  bio-optical  products 
on  a  daily  basis  and  provide  a  new  capability  for  monitoring  coastal  processes.  A  visual  comparison  of  our  new 
composite  product  against  our  previous  operational  composite,  the  7-day  latest  pixel  composite,  shows  that  the 
new  process  yields  a  cleaner,  more  representative  picture  of  the  oceanographic  fields. 
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Figure  7.  The  merge  and  interpolation  of  the  7-day  LPC  and  the  previous  three  steps. 
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Figure  8.  The  final  step  of  using  the  nearest  neighbor  fill  results  in  a  completely  data 


filled  chlorophyll  image 
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Cloud  contaminated  pixels  are  still  an  issue  and  produce  some  compositing  artifacts.  Additionally,  selecting 
an  overly  large  search  radius  can  produce  strange  looking  results  when  the  triangular  interpolation  is  performed 
over  large  distances. 
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